#'###########################################
#'                                          #
#'  ------------------------------------    #
#'  |     COIN     __                  |    #
#'  |            >(' )                 |    #
#'  |              )/                  |    #
#'  |             /  `----/            |    #
#'  |             \  ~=- /             |    #
#'  |   ~^~^~^~^~^~^~^~^~^~^~^~^~^~^   |    #
#'  ------------------------------------    #
#'                                          #
#'###########################################

rm(list=ls())

source("C:/UQAR/Recherche/Maitrise/Donnees/R/Utils/loadDependencies.R")


STATIC_WORKING_DIR <- "C:/UQAR/Recherche/Maitrise/Donnees/Donnees/Bylot/Fixe"


#'#
#'# Recuperation des donnees statiques et comparaison avec le gps
#'#

static.errors.filtered <- loadData(STATIC_WORKING_DIR)
static.errors.unfiltered <- loadData(STATIC_WORKING_DIR,
				locsclasses=c(3, 2, 1, 0, "A", "B", "Z"))


#'#################################################################################################
#'#                                                                                              ##
#'#                                 Proportion de localisations                                  ##
#'#                                                                                              ##
#'#################################################################################################



#'#
#'# Compositions
#'# 


static.prop.test <- prop.table.byLC(static.errors.filtered)
static.prop.test$session <- rep(1:5, each=3)
replaceall.na.dt(static.prop.test, NaN)

fact <- static.prop.test[, list(visibility=as.factor(visibility), id=as.factor(id), session=as.factor(session))]
static.prop.test$visibility <- NULL
static.prop.test$id <- NULL
static.prop.test$session <- NULL

imp.comp <- impCoda(static.prop.test, k=3, method="lm")$xImp
static.comp <- acomp(imp.comp)
#static.comp <- acomp(static.prop.test, BDL=NaN)
summary(static.comp)

(mylm <- lm(alr(static.comp) ~ visibility, data=fact))
summary(manova(mylm), test="Wilks")



datat <- cbind(clr(static.comp), as.data.frame(fact))
colnames(datat)[1:5] <- c("LC3", "LC2", "LC1", "LC0", "LCA")
dmelt <- melt(datat, id.vars=c("visibility", "id", "session"))


lmLc3 <- lm(LC3 ~ visibility + id, data=datat)
anova(lmLc3)
lmLc2 <- lmer(LC2 ~ visibility + (1|id), data=datat)
anova(lmLc2)
lmLc1 <- lmer(LC1 ~ visibility + (1|id), data=datat)
anova(lmLc1)
lmLc0 <- lmer(LC0 ~ visibility + (1|id), data=datat)
anova(lmLc0)
lmLca <- lmer(LCA ~ visibility + (1|id), data=datat)
anova(lmLca)

mlm <- lmer(value ~ -1 + variable + visibility + (variable - 1|id), data=dmelt)
summary(mlm)
anova(mlm)
difflsmeans(mlm)


mm <- MCMCglmm(cbind(LC3, LC2, LC1, LC0, LCA)~ visibility, random=~id, data=datat,
    rcov=~us(trait):units, family=c("categorical","categorical"))


#'#################################################################################################
#'#                                                                                              ##
#'#                                 Evaluation de l'erreur                                       ##
#'#                                                                                              ##
#'#################################################################################################


r <- static.errors.filtered[, list(mean=mean(distance), q=quantile(distance, .68)), by=list(visibility)]
(r$q[1] - r$q[3]) / r$q[1] * 100
(r$q[2] - r$q[3]) / r$q[2] * 100
(r$mean[1] - r$mean[3]) / r$mean[1] * 100
(r$mean[2] - r$mean[3]) / r$mean[2] * 100


#'#
#'# Effet de la visibilite : Modele mixte
#'#

aovdata <- static.errors.filtered[, list(visibility, LC, id, distance, diffLong, diffLat)]
aovdata[, visibility := as.factor(visibility)]
aovdata[, id := as.factor(id)]

lmm <- lmer(log(distance) ~ visibility * LC + (1|id), data = aovdata)
summary(lmm)
difflsmeans(lmm)

x11()
plot(lsmeans(lmm))

plot(difflsmeans(lmm))

anova(lmm)


plotLMER.fnc(lmm, pred="LC", intr=list("visibility", c("1", "2", "3"), "beg", 
				list(c("black", "red", "yellow"), rep(1,3))))


randomEffects <- ranef(lmm, postVar = TRUE)
qqmath(randomEffects)
dotplot(randomEffects)
rsquared.lme(list(lmm))




#'#
#'# Distributions cumulatives
#'#
													
cumuldist <- list()
for (i in locsclasses) {
	cumuldist[i] <- list(ecdf(static.errors.filtered[i][, distance]))
}

localisations.props <- do.call(rbind, lapply(cumuldist, prop.between, props=c(250, 500, 1500)))
colnames(localisations.props) <- c("0-250", "250-500", "500-1500", ">1500")


colors <- rainbow(6)
names(colors) <- locsclasses
#plot(names(res), res, log="x")
plot(cumuldist[["3"]], log="x", xlim= c(10, 100000), do.points=F, verticals=T, col=colors[1], lwd=2)
for (i in 2: length(cumuldist)) {
	plot(cumuldist[[i]], add=TRUE, do.points=F, verticals=T, col=colors[i], lwd=2)
}
abline(v=250)
abline(v=500)
abline(v=1500)
legend("bottomright", col=colors, legend=names(cumuldist))

barplot(localisations.props, beside=TRUE)
barplot(t(localisations.props), beside=TRUE)


